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ABSTRACT:  Angles-only  (or  bearings-only)  navigation  involves  determining  position,  velocity,  or  orientation  in¬ 
formation  for  an  observer  using  the  apparent  directions  or  motions  of  objects  at  finite  distances.  Angles-only  navi¬ 
gation  covers  a  broad  range  of  applications,  and  is  generally  implemented  through  a  Kalman  filter  that  uses 
imaging  and  other  information  to  differentially  adjust  the  values  of  the  navigation  parameters  (state  vector)  at 
each  incremental  step  of  the  observer’s  computed  motion. 

This  paper  presents  an  algorithm  for  angles-only  navigation  that  is  a  closed-form  solution  for  both  position  and 
velocity  that  does  not  require  any  prior  estimate  of  the  observer’s  position  or  motion.  It  is  least-squares-based  trian¬ 
gulation  generalized  to  a  moving  observer,  involving  only  angular  observations  of  objects  with  known  coordinates. 

The  algorithm  can  be  applied  to  any  situation  where  foreground  objects  are  observed  against  background 
objects,  so  long  as  both  foreground  and  background  objects  are  identifiable  and  have  known  coordinates.  It  can 
be  used  for  ship  piloting  using  images  of  shore  objects  at  various  distances.  It  also  has  possible  applications  in 
robotics;  for  example,  for  determining  the  position  and  velocity  of  mobile  landers  on  other  solar  system  bodies,  or 
for  automated  farming.  A  specific  proposed  application  for  open-sea  ship  navigation  would  use  the  angular  posi¬ 
tions  of  Earth  satellites  observed  optically  against  a  star  background.  Such  a  system  could  provide  a  supplement 
to  ordinary  GPS  navigation,  as  well  as  supplying  a  precise  absolute  attitude  reference. 


INTRODUCTION 

This  paper  outlines  the  geometry  and  correspond¬ 
ing  mathematics  of  a  particular  type  of  angles-only 
navigation.  In  angles-only  navigation  (also  known  as 
bearings-only  navigation)  position,  velocity,  or  orien¬ 
tation  information  for  an  observer  is  passively  ob¬ 
tained  from  measurements  of  the  apparent  angles, 
or  angular  rates,  of  objects  at  finite  (but  generally 
unknown)  distances.  Despite  its  name,  angles-only 
navigation  is  often  used  to  augment  other  means  of 
navigation,  such  as  dead  reckoning,  inertial,  or  GPS. 
As  a  simple  example,  a  ship’s  navigator  can  deter¬ 
mine  a  line  of  position  from  the  measured  bearing  of 
an  identifiable  navigation  aid  or  landmark.  Two  such 
lines  of  position  will  cross  at  the  ship’s  position.  In 
more  sophisticated  forms,  it  has  been  applied  to 
spacecraft  maneuvering  [1,  2],  aircraft  navigation 
[3-6],  and  position  determination  for  mobile  robots 
[7-9].  Advances  in  electronic  imaging  and  real-time 
image  analysis  capabilities  over  the  last  few  decades 
have  considerably  expanded  the  scope  of  uses,  and 
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the  literature  base  has  blossomed,  with  significant 
contributions  from  fields  such  as  optics,  computer 
vision,  robotics,  artificial  intelligence,  and  even  ani¬ 
mal  navigation.  The  topic  is  probably  too  broad  for 
a  review  even  to  be  possible,  and  the  references 
cited  above  are  just  a  tiny  sample  of  the  papers 
published. 

In  the  most  common  implementations  of  angles- 
only  navigation,  measurements  from  a  scene  re¬ 
corded  by  an  imaging  system  serve  as  input,  along 
with  data  from  other  sensors,  to  a  navigational 
Kalman  filter  that  continually  updates  the  observ¬ 
er’s  state  vector  (e.g,  position,  velocity,  and  atti¬ 
tude);  see,  for  example,  [10].  For  this  purpose,  the 
sensitivity  of  the  scene  elements  and  other  sensor 
data  to  a  change  of  state  are  linearized  about  an 
estimated  state,  which  is  a  computational  projec¬ 
tion  based  on  previous  data.  The  differences  be¬ 
tween  the  measurements  and  their  expected  values, 
which  are  assumed  small,  provide  information  to 
correct  the  estimated  state,  and  the  cycle  repeats. 
A  somewhat  different  kind  of  predictive  stepwise 
filter  for  angles-only  measurements,  yielding  posi¬ 
tion,  velocity,  attitude,  and  rotation,  was  published 
in  [11], 

In  this  paper,  however,  we  consider  how  to  deter¬ 
mine,  ab  initio,  both  the  position  and  velocity 
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vectors  of  an  observer  (e.g.,  an  imaging  system) 
using  a  sequence  of  angular  measurements  that  are 
distributed  over  a  period  of  time.  The  measure¬ 
ments  are  of  the  apparent  directions  of  identifiable 
objects  with  known  coordinates.  The  scheme  is 
“absolute”  in  the  sense  that  the  angular  measure¬ 
ments,  rather  than  being  relative  to  some  unknown 
directions,  are  expressed  in  the  reference  system 
that  is  used  for  the  object  coordinates  (and  the  navi¬ 
gation  solution).  This  considerably  simplifies  the 
problem  by  obviating  the  need  for  a  simultaneous 
attitude  solution  or  any  kind  of  object-space/image- 
space  mapping.  Although  at  first  thought  this  situa¬ 
tion  may  seem  atypical,  the  development  remains 
pertinent  to  many  common  applications.  For  exam¬ 
ple,  it  is  relevant  to  any  system  that  captures 
scenes  in  which  foreground  features  appear  against 
background  features,  and  geodetic  coordinates  can 
be  obtained  for  both  near  and  far  objects. 

In  contrast  to  the  Kalman  filter  approach,  the 
algorithm  presented  here  does  not  require  any  pre¬ 
vious  estimate  of  position  or  motion,  and  is  of  closed 
form,  not  stepwise  or  iterative.  It  is  a  type  of  3-D 
triangulation  applied  to  a  moving  observer,  with 
angular  observations  taken  at  various  positions 
along  the  observer’s  track.  The  observations  are 
assumed  to  be  uncorrelated  and  to  have  normally 
distributed  random  errors  but  no  significant  sys¬ 
tematic  errors.  The  solution  minimizes  the  effects 
of  errors  in  both  the  observations  and  the  assumed 
object  coordinates  in  a  least-squares  sense.  Since 
the  algorithm  requires  only  information  related  to 
the  angular  observations,  it  may  be  useful  for 
startup  situations  or  any  other  circumstances  in 
which  current  position  and  velocity  data  are  unreli¬ 
able  or  not  available.  The  scheme  requires  only 
that  the  objects  observed  can  be  identified  and  their 
coordinates  retrieved.  The  algorithm  can  also  pro¬ 
vide  a  check  on  stepwise  navigation  filters.  The  de¬ 
velopment  is  based  on  a  straightforward  geometric 
construction  and  the  solutions  are  robust  for  rea¬ 
sonable  sets  of  observations. 

•  What  is  presented  here  is  not,  of  course,  the  first 
closed-form  solution  to  an  angles-only  navigation 
problem.  In  fact,  the  solution  for  the  fixed-observer 
case  was  published  (in  a  different  kind  of  notation) 
in  an  appendix  to  the  classical  text  Geodesy  by 
Bomford  [12].  Two  of  the  robotics  papers  mentioned 
above  [7,  9]  present  closed-form  solutions  for  a  fixed 
observer  using  relative  bearing  measurements  in  a 
2-D  environment.  The  main  contribution  of  this  pa¬ 
per  is  in  presenting  a  closed-form  solution  for  both 
position  and  velocity  in  a  3-D  environment.  The  de¬ 
velopment  includes  a  correction  term  for  the  curva¬ 
ture  of  the  Earth,  so  that  observations  can  be  col¬ 
lected  over  extended  tracks. 

The  paper  is  organized  by  sections  as  follows. 
First,  some  possible  applications  of  the  method  are 


described.  The  following  section  explains  the  basic 
concepts  of  the  observations  and  their  representa¬ 
tion  as  vectors.  Next,  two  navigation  solutions  are 
presented,  one  for  a  fixed  observer  and  the  other  for 
a  moving  observer.  The  term  in  the  solution  that 
corrects  for  the  curvature  of  the  Earth  is  then 
briefly  discussed.  The  next  two  sections  consider 
the  propagation  of  random  error  and  the  effects  of 
an  unmodeled  acceleration  on  the  solution.  The  spe¬ 
cific  possibility  of  using  Earth  satellites  observed 
optically  (or  in  the  near-infrared)  against  back¬ 
ground  stars  as  the  source  of  observations  for  this 
method  is  described  in  some  detail.  The  results  of  a 
numerical  simulation  of  such  an  observing  system 
are  presented  and,  finally,  the  salient  points  of  the 
paper  are  summarized. 


APPLICATIONS 

The  method  described  in  this  paper  was  developed 
for  a  proposed  shipboard  automated  celestial  observ¬ 
ing  system  that  would  augment  GPS  with  absolute 
orientation  information  and  serve  as  a  standalone 
positioning  system  in  case  of  GPS  denial.  The  cho¬ 
sen  observing  mode  involved  artificial  Earth  satel¬ 
lites  observed  optically  or  in  a  near-infrared  band 
against  background  stars.  The  algorithm  described 
here  was  initially  developed  simply  to  provide  quick 
estimates  of  the  likely  errors  of  such  a  system, 
under  various  conditions,  even  though  the  eventual 
navigation  solution  was  to  be  obtained  from  a  Kal¬ 
man  filter  involving  multiple  sensor  inputs.  How¬ 
ever,  the  method  has  value  in  itself  by  providing  a 
standalone  navigation  solution  from  the  satellite 
observations  in  the  absence  of  any  other  informa¬ 
tion.  The  obvious  next  question  was  whether  the 
algorithm  might  be  applied  in  other  contexts. 

As  it  turns  out,  the  mathematics  provided  here 
could  be  applied  to  any  situation  in  which  the  direc¬ 
tions  to  identifiable  objects  can  be  measured  with 
respect  to  more  distant  objects,  so  long  as  the  coor¬ 
dinates  of  both  foreground  and  background  objects 
are  known.  When  a  foreground  and  a  background 
object  appear  to  line  up  from  the  point  of  view  of 
the  observer,  the  observed  direction  vector  is  simply 
the  normalized  difference  of  the  position  vectors 
of  the  two  objects.  Since  the  direction  vector  is  com¬ 
puted  from  absolute  coordinates  —  that  is,  coordi¬ 
nates  relative  to  a  well-defined  reference  system, 
such  as  WGS-84  —  it,  too,  is  absolute.  Using  that 
type  of  observation,  the  algorithm  therefore  has 
wide  application  to  any  kind  of  vehicle  that  has 
an  imaging  system,  even  if  not  of  the  highest  qual¬ 
ity.  The  accuracy  depends  only  on  the  resolution 
of  the  image,  not  on  any  external  angular  calibra¬ 
tion  or  the  transformation  of  real  space  to  image 
space. 
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Fig.  1-Images  of  aligned  foreground  and  background  objects  viewed  from  the  Patapsco  River  out¬ 
side  Baltimore.  Each  image  yields  a  direction  vector  (line  of  position)  in  the  WGS-S4  system,  based 
on  the  coordinates  of  the  two  objects. 


For  objects  seen  near  the  horizontal  plane,  the 
apparent  alignment  of  objects  in  the  vertical  direc¬ 
tion  suffices  if  the  height  of  the  observer  is  not  of  in¬ 
terest  (i.e.,  if  the  3-D  problem  effectively  collapses 
to  2-D).  For  example,  as  seen  from  a  ship,  a  buoy 
may  appear  to  pass  underneath  a  distant  water 
tower.  A  series  of  time-tagged  images  of  aligned 
objects,  for  which  the  latitudes  and  longitudes  are 
known,  would  allow  for  a  solution  for  the  surface 
track  of  the  vehicle.  The  method  could  be  applied  to 
ordinary  ship  piloting,  which  was  demonstrated  by 
the  author  using  images  of  shore  objects  taken  dur¬ 
ing  a  small-boat  trip  up  the  Chesapeake  Bay  and 
into  Baltimore  Harbor  (see  Figure  1  for  examples). 
The  images  were  taken  with  an  inexpensive  hand¬ 
held  camera  and  timed  only  to  the  nearest  minute 
using  the  camera’s  internal  clock.  For  each  image, 
coordinates  of  the  foreground  and  background 
objects  were  subsequently  obtained  from  Google 
Earth  or  the  USCG  Light  List  [13],  and  the  coordi¬ 
nates  converted  to  geocentric  position  vectors.  Using 
two  vectors  derived  from  each  such  observation  (see 
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next  section),  a  solution  for  a  portion  of  the  boat’s 
track  was  then  computed  using  the  method  de¬ 
scribed  in  this  paper.  Despite  the  primitive  nature 
of  the  experiment,  the  solution  was  close  to  that 
obtained  from  a  straight-line  fit  of  the  recorded  GPS 
positions  of  the  boat,  and  the  residual  errors  in  both 
cases  (748  and  629  m  RMS,  respectively)  were 
largely  the  result  of  the  boat’s  deviations  from  the 
modeled  tracks  —  an  important  consideration  in  the 
practical  applicability  of  the  method  that  will  be  dis¬ 
cussed  later.  This  rudimentary  exercise  showed, 
however,  that  the  method  might  be  useful  for  pilot¬ 
ing  in  areas  that  are  well  mapped,  or  for  which  good 
satellite  imagery  is  available,  but  lack  reliable  navi¬ 
gational  aids.  Practical  use  of  the  method  would 
require  an  automated  system  of  object  identification 
that  would  permit  the  retrieval  of  their  coordinates 
from  a  suitable  database. 

The  method  would  also  apparently  work  well  for 
navigating  robotic  landers  across  the  surfaces  of 
other  solar  system  bodies  using  the  on-board  imaging 
system.  All  that  is  needed  is  a  database  of  identifiable 
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Fig.  2-Geometry  of  a  single  observation.  Both  the  observed  direction  of  the  object  and  the  object’s 
coordinates  are  assumed  to  have  some  error. 


terrain  features  and  their  coordinates.  In  fact,  it 
could  be  applied  to  interplanetary  space  navigation 
(although  not  at  high  accuracy)  if  the  apparent 
directions  of  several  relatively  nearby  solar  system 
objects  were  observed  against  the  star  background. 
Navigation  of  automated  farming  systems  is  another 
possibility.  For  example,  a  system  of  marked  poles 
at  a  field’s  edge,  at  two  different  perimeters,  might 
be  imaged  by  a  rotating  camera  on  top  of  a  robotic 
tractor. 


OBSERVATIONS,  VECTORS,  COORDINATE 
SYSTEMS 

This  paper  uses  the  convention  that  vectors  of  ar¬ 
bitrary  length  are  written  as  boldface  upper-case 
letters  and  unit  vectors  are  written  as  boldface 
lower-case  letters.  For  example,  z  would  be  the  unit 
vector  in  the  direction  of  Z. 

The  algorithm  in  this  paper  is  based  on  observa¬ 
tions  of  the  directions  of  identifiable  objects,  with 
known  coordinates,  from  the  point  of  view  of  an  ob¬ 
server  whose  own  coordinates  are  to  be  determined. 
For  each  object  observed,  then,  two  kinds  of  infor¬ 
mation  are  required:  the  predetermined  coordinates 
of  the  object,  represented  by  the  position  vector  P; 
and  the  observation  itself,  represented  by  the  direc¬ 
tion  (unit)  vector  d.  In  the  absence  of  errors,  the 


observer  must  be  somewhere  on  a  line  of  position 
(LOP)  in  3-D  space  given  by  the  equation  X  =  P  + 
rd,  where  X  is  the  position  of  an  arbitrary  point 
along  the  line  and  r  is  a  scalar  that  can  take  on  any 
real  value.  The  components  of  the  vectors  X  and  P 
and  the  scalar  r  have  units  of  length,  while  d  is 
dimensionless.  We  assume  that  X,  P,  and  d  may  be 
functions  of  time;  for  a  moving  target,  the  time  se¬ 
ries  of  vectors  P it)  is  referred  to  as  its  ephemeris. 
See  Figure  2. 

We  require  that  P  and  d  are  given  in,  or  reduced 
to,  a  common  coordinate  system.  For  some  kinds  of 
observations,  this  will  come  about  naturally.  For 
example,  if  a  target  object  is  observed  against  a 
background  object,  and  both  have  coordinates  in 
the  same  database,  then  the  direction  vector  is  sim¬ 
ply  the  difference  between  the  known  position  vec¬ 
tors  of  the  background  object  and  the  target,  nor¬ 
malized  to  unit  length.  The  direction  vector  is 
thereby  expressed  in  the  coordinate  system  used  for 
the  positions  of  all  the  landmarks.  A  more  compli¬ 
cated  case  is  that  of  artificial  Earth  satellites 
imaged  against  the  star  background,  in  which  the 
observations  and  object  coordinates  are  naturally 
expressed  in  different  Icinds  of  coordinates  and 
must  be  reduced  to  a  common  system.  This  is  dis¬ 
cussed  more  fully  in  the  section  on  using  satellites 
as  target  objects. 
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The  navigation  solution  vectors  will  be  expressed 
in  the  coordinate  system  used  for  the  object  coordi¬ 
nates  and  observations.  Effectively,  the  set  of 
assumed  coordinates  for  all  the  objects  observed 
define  the  reference  system  for  the  solution.  Since 
the  solution  is  based  on  geometry,  not  dynamics, 
there  is  no  requirement  that  the  reference  system 
be  inertial.  In  the  most  common  case,  in  fact,  the 
reference  system  will  be  geodetic  (Earth-fixed)  and 
therefore  not  inertial. 


NAVIGATION  SOLUTION 

We  will  start  with  an  observer  that  is  fixed  with 
respect  to  the  objects  he  is  observing,  at  location  X. 
Given  an  ensemble  of  n  observations  d,  of  objects  at 
known  positions  P;,  respectively,  the  least-squares 
estimate  for  X  is  given  by 


n~[  dfj 

[djjdjj] 

-[d;id;J  ^ 

/  xi 

[d/i  dj2] 

»  - 

[d;2djj] 

X2 

-[c MU 

\x3 

/{P£l-(d/.Pi)dil]\ 

=  [P,-t  -  (df  ■  Pi)d*a]  (1) 

V  [Pia  -  (d;  •  P;)di3]  ) 

where  P;  =  (P^,  P,2,  Pq),  d,  =  (dfi,  d;2,  dis),  X  =  (x1;  x2, 
X3),  and  the  square  brackets  indicate  a  summation 

n 

over  all  n  observations,  i.e.,  [. . .]  represents  J2  ■  ■  ■■ 

i= 1 

This  is  a  system  of  three  scalar  equations  in  three 
unknowns,  xi,  x2,  and  X3,  the  components  of  the 
observer’s  position  vector.  These  three  equations 
are  equivalent  to  Equation  (C.29)  in  [12].  We  have 
minimized  the  sum  V  -Yhi  where  each  repre¬ 
sents  the  distance  of  line  of  position  i,  defined  by 
the  observation  i,  from  X,  the  computed  position  of 
the  observer.  We  imagine  all  of  the  LOPs  converg¬ 
ing  in  a  small  volume  of  3-D  space;  X  is  in  the  cen¬ 
ter  of  this  volume,  with  “center”  precisely  defined 
by  the  least-squares  criterion. 

There  is  a  brute-force  solution  to  Equation  (1) 
from  substitution.  If  we  re-cast  (1)  as 

(A  B  C\  fx1\  f Qi\ 

B  D  E  x2  =  Q2 

\C  E  Fj  VW  \Qj 

with  A  =  n  -  [d“],  B  =  “[d(,dj,  etc.,  then  the  solu¬ 
tion  is 


{CD-BE){CQ-Sl  -AQ3)  -  (. BC  - AE){CQ2-BQ3 ) 
{CD  - BE){C 2  -AF)-  {BC -AE){CE- BF) 
CQ2-BQ3~{CE-BF)x3 
CD -BE 

Qi  Bx  2  Cx  3  f  o\ 


Because  this  development  applies  only  to  a  fixed 
observer,  or  n  simultaneous  observations,  it  is  of 
limited  use  for  navigation.  For  simplicity,  we  have 
also  not  considered  weighting  the  observations, 
although  providing  for  that  is  straightforward,  and 
is  described  at  the  end  of  this  section. 

The  more  useful  case  is  that  of  a  moving  observer 
and  non-simultaneous  observations.  Let  us  repre¬ 
sent  the  observer’s  trajectory  by  X(t)  —  X(£0)  + 
V{t0)t  4-  /h)x0,  where  X(t.0)  and  V(t0)  are  the  observ¬ 
er’s  position  and  velocity,  respectively,  at  time  t0.  In 
the  third  term,  x0  is  the  unit  vector  in  the  direction 
X(t0)  and  fit)  is  a  scalar  function  with  units  of  dis¬ 
tance.  The  time  t  is  measured  from  the  time  origin 
to,  which  can  be  chosen  for  convenience;  it  is  an  ar¬ 
bitrary  time  that  is  within  or  not  far  outside  the 
span  of  observation  times  (ii  to  tn)  but  it  does  not 
necessarily  correspond  to  the  time  of  a  specific  ob¬ 
servation.  The  observations  occur  at  discrete  times 
ti  (measured  from  i0),  and  for  those  times  the  trajec¬ 
tory  can  be  expressed  as  Xj  =  X0  +  VQ^  4-  /)x0, 
using  the  shorthand  X*  =  X(tf,  Xo  =  X(t 0),  V0  = 
V(t0),  and  ft  =  ftf. 

The  third  term  represents  any  curvature  in  the 
observer’s  path  in  the  direction  x0,  which,  if  a  geo¬ 
centric  coordinate  system  is  used  and  fit)  <  0,  is  to¬ 
ward  the  center  of  the  Earth.  Thus,  the  third  term 
could  represent  the  gravitational  acceleration  of  an 
object  in  Earth  orbit  or,  for  an  observer  traveling  on 
or  near  Earth’s  surface,  the  local  curvature  of  the 
geoid.  If  the  third  term  is  written  as  fjXfJXo,  where 
X0=  !  Xq  I ,  then  X;  =  X0(l  +|^)  +  V0£;.  The  curva¬ 
ture  term  is  assumed  small  compared  to  the  other 
terms,  i.e.,  that  fi  <C  Xo  and  f  <  1  Vo  I  Ai,  where 
At  =  I  tn  -  tj  I  is  the  span  of  time  covered  by  the 
observations.  It  is  also  assumed  that  f/Xo  (which  is 
small)  can  be  considered  known  to  sufficient  accu¬ 
racy.  The  calculation  of  ft  can  be  done  in  a  number 
of  ways  and  is  not  essential  to  the  method;  see  next 
section. 

The  equation  below  represents  a  least-squares  so¬ 
lution  to  the  navigation  problem;  specifically,  it 
minimizes  the  sum  V  —  bf  where  each  <5;  repre¬ 
sents  the  distance  of  line  of  position  i,  defined  by 
the  observation  taken  at  time  ti,  from  X/,  the  com¬ 
puted  position  of  the  observer  at  the  same  instant. 
The  geometric  interpretation  of  the  solution  is  simi¬ 
lar  to  that  for  the  fixed-observer  case,  but  less  intui¬ 
tive.  The  LOPs  do  not  converge  to  define  a  small 
volume  of  space;  rather,  they  converge  around  the 
observer’s  computed  path  in  such  a  way  that  at  the 
time  of  each  observation,  the  observer  is  as  close  as 
possible  to  the  LOP  given  the  simple  model  we  are 
adopting  for  his  motion  and  the  fact  that  the  close¬ 
ness  criterion  (square  of  the  distance)  is  assessed  in 
the  aggregate  for  all  the  LOPs. 

The  solution  algorithm,  derived  elsewhere  [14] 
is: 
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where,  again,  the  square  brackets  indicate  a  sum¬ 
mation  over  all  n  observations,  and  ft  — 

The  column  vector  on  the  left  side  represents  the 
unknown  navigation  state,  U,  at  time  t0:  U  =  (Xo, 
Vo)  =  (xi,  X2,  x3,  Vi,  v2,  v3).  The  remaining  notation 
is  the  same  as  for  Equation  (1).  Note  that  the  time 
A  of  each  observation,  measured  from  *o>  appears 
explicitly  in  some  of  the  terms.  The  time  *o  can  be 
chosen  to  be  at  the  end  of  the  series  of  observations, 
so  that  Equation  (3)  can  provide  a  near  real-time 
estimate  of  position  and  velocity. 

Equation  (3)  is  of  the  form  AU  =  Q,  where  U  and 
Q  are  column  6-vectors  and  A  is  a  6x6  matrix.  The 
solution  is  U  =  A_1Q,  where  A-1  is  the  inverse  of  A 
and  represents  the  unsealed  covariance  matrix  of 
the  solution. 

For  straight-line  motion  (no  curvature  term), 
ft= 1  for  all  i.  Also,  if  the  observer  is  stationary, 
then  tt  can  be  considered  to  be  0  for  all  i  since  time 
is  measured  from  when  the  observer  was  at  Xq  (the 
problem  then  is  three-dimensional  rather  than  six¬ 
dimensional).  In  that  case,  (3)  reduces  to  (1).  If  the 
observer  is  moving  but  the  velocity  vector  is  known, 
the  position  vector  can  be  obtained  from  the  first 
three  rows  on  the  left  side  of  (3)  (actually,  any  three 
rows)  if  the  terms  involving  v1;  v2,  and  v3  are 
moved  to  the  right  side. 

If  the  observations  have  different  uncertainties, 
then  the  algorithm  should  minimize  the  weighted 
sum,  T>w  =  Yjiiwi^i)2>  where  Wi  is  the  dimension¬ 
less  weight  of  observation  i.  The  weight,  wt  =  c/cy 
is  the  ratio  of  the  average  uncertainty  of  all  the 
observations,  a,  to  the  uncertainty  of  the  particular 
observation,  ay  (Often  it  is  desirable  that  J2  wf  =  n, 
so  that  1/a2  =  (l/crf).)  The  uncertainty  of  an  ob¬ 

servation  is  usually  dominated  by  the  angular 
measures  that  define  the  vector  d,  for  the  observa¬ 
tion  (see  section  below  on  random  errors).  Including 


observational  weights  is  accomplished  simply  by 
including  the  extra  factor  wf  in  each  of  the  sums 


in  (3). 

lijif 


CURVATURE  TERM 


The  third  term  in  our  motion  model  X(*)  —  X(f0)  4- 
V0f  +  flt)x0  describes  the  curvature  of  the  path  in 
the  direction  x0,  which,  for  a  geocentric  coordinate 
system  and  fit )  <  0,  will  be  toward  the  center  of  the 
Earth.  The  term’s  value  at  the  time  of  observation  i 

is  fi  —  fitf))  the  term  can  then  be  written  (^JQXq  and 

we  have  assumed  that  ft/X0  is  a  known  (dimension¬ 
less)  quantity.  For  short  tracks  on  the  surface  of  the 
Earth,  ft/X 0  is  small.  For  example,  even  for  a  100- 
km  track,  ft  =  -0.8  km,  so  fJX '0  «  10-4. 

In  many  cases,  the  curvature  term  may  be 
unnecessary.  It  would  be  important  for  navigation 
applications  on  or  near  the  surface  of  the  Earth  in 
which  observations  are  collected  over  a  track  that 
may  extend  to  some  tens  of  kilometers  or  more  and 
navigational  accuracies  of  better  than  100  m  are 
expected.  Generally  these  would  be  aircraft  applica¬ 
tions.  It  is  used  in  place  of  an  acceleration  term 
in  the  solution,  which  would  require  three  more 
unknowns  (and  at  least  two  more  observations)  and 
which  would,  in  many  cases,  be  poorly  determined 
because  of  its  smallness  compared  to  observational 
error. 

For  most  purposes  the  magnitude  of  the  term  can 
be  well  represented  by  the  parabolic  approximation 


1  {vtjf 

2  R 


or 


fi  _  1  {vtjf  =  _  1  / vtf\ 

Xo  2  RX 0  2  / 
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where  R  is  the  radius  of  curvature  of  the  path,  v  — 

I  Vo  I  is  the  speed  of  motion,  and  oft  <g;  R.  The  length 
vti  is  the  distance  traveled  from  the  reference  point 
(where  t  —0  and  X  =  X0)  to  the  position  of  observation 
i,  taken  at  time  A-  For  an  observer  in  a  circular  Earth 
orbit,  R  =  Xq  and  v2  =  GM/X0,  where  GM  is  the  geo¬ 
centric  gravitational  constant,  but  we  will  not  con¬ 
sider  the  orbital  case  further  here.  For  an  observer 
on  or  near  the  surface  of  the  Earth,  a  great-circle 
course  is  implied.  For  such  an  observer,  if  we  con¬ 
sider  the  Earth  to  be  a  sphere,  R  —  Xq  —  a  +  h,  where 
a  is  the  radius  of  the  Earth  and  h  is  the  height  above 
sea  level.  Note  that  some  pre-solution  estimate  of 
speed  is  necessary  for  the  evaluation  of  the  term. 

On  the  real,  oblate  Earth,  things  are  a  bit  more 
complicated,  because  R  A  Xo,  that  is,  the  local  ra¬ 
dius  of  curvature  is  not  the  same  as  the  local  geo¬ 
centric  distance,  and  both  vary  from  place  to  place. 
Formulas  for  radius  of  curvature  and  distance  to 
the  geocenter  are  given  in  any  elementary  geodesy 
text,  for  example,  [12,  15,  16].  However,  in  most 
cases,  the  spherical-Earth  approximation  will  work 
sufficiently  well.  For  high-accuracy  applications,  or 
those  where  the  observation  track  is  extended,  only 
very  crude  estimates  of  the  observer’s  position, 
direction,  and  speed  are  needed  to  evaluate  the 
term.  The  curvature  term  is  discussed  more  fully  in 
[14],  which  provides  estimates  of  its  sensitivity  to 
assumed  location  and  speed. 


PROPAGATION  OF  RANDOM  ERRORS 

As  an  application  of  least-squares,  the  algorithm 
represented  by  Equation  (3)  assumes  that  the 
observations  are  uncorrelated  and  that  they  have 
normally  distributed  random  errors.  Systematic 
errors  are  assumed  to  be  insignificant.  The  effects 
of  some  common  systematic  errors  are  treated  in 
the  next  section.  This  section  discusses  the  origin 
and  propagation  of  random  errors  of  measurement. 

The  algorithm  given  above  differs  from  most 
least-squares  applications  in  two  important  ways. 
First,  the  quantity  that  is  minimized  in  a  sum-of- 
squares  sense  is  a  euclidean  distance  in  3-D  space 
that  is  related  to,  but  is  not  itself,  a  measured 
quantity  for  each  observation.  These  distances  (ft) 
will  have  a  statistical  scatter  that  reflects  not  just 
the  errors  of  measurement  but  also  the  errors 
in  the  predetermined  coordinates  of  the  objects 
observed.  Second,  the  method  does  not  rely  on  a  lin¬ 
earization  around  an  approximately  known  set  of 
parameters.  The  solved-for  parameters  are  the  posi¬ 
tion  and  velocity  vector  components,  not  corrections 
to  the  components  of  assumed  vectors.  Yet,  despite 
the  fact  that  no  conditional  (observation)  equations 
have  been  defined,  most  aspects  of  least-squares 
f  analysis  still  apply. 


For  example,  as  previously  stated,  the  inverse 
of  the  6X6  matrix  in  Equation  (3),  A-1,  is  the 
unsealed  covariance  matrix  of  the  solution.  We  can 
use  it  in  the  conventional  way  to  obtain  the  formal 
uncertainties  of  the  six  solved-for  parameters,  (x1; 
x2,  x3,  Vj,  v2,  v3),  and  the  parameter  correlation 
matrix: 


where  a?  is  the  formal  variance  of  parameter  i  (i= 1 
to  6),  fffj  is  the  formal  covariance  of  parameters  i 
and  j,  and  c ij  is  the  correlation  (—1  to  +1)  between 
parameters  i  and  j.  We  could,  then,  use  the  varian¬ 
ces  and  covariances  of  Xi,  x2,  and  x3  to  form  an 
error  ellipsoid  in  3-D  space  that  represents  the 
uncertainties  in  the  solution  for  X0,  or  propagate 
the  formal  errors  to  any  other  position  on  the 
observer’s  track  to  see  how  the  error  ellipsoid 
changes  with  time. 

The  factor  Vl{2n  -  6)  represents  the  variance  of 
the  least-squares  fit,  i.e.,  a  measure  of  the  scatter 
in  the  post-solution  residuals.  We  have  V  =  ]CA  > 
where  each  ft  is  the  distance  between  the  LOP  of 
observation  i  and  the  computed  position  of  the  ob¬ 
server  at  the  same  time.  The  quantity  2 n  -  6  repre¬ 
sents  the  number  of  degrees  of  freedom  in  the  solu¬ 
tion  and  reflects  the  fact  that  each  of  the  n  observa¬ 
tions  is  two-dimensional,  i.e.,  it  consists  of  two 
independent  angular  measurements.  Each  di  in  the 
sum  can  be  calculated  using 

ft  =  idjX(Pi  -X)|  (6) 

which  is  the  minimum  distance  from  point  X;  to  its 
associated  LOP,  P,  +  r<ft.  The  symbols  have  the 
same  meanings  as  in  Equations  (1)  and  (3),  and 

X£  =  X0  (l  +  J-j  +  V0A,  with  Xo  and  V0  taken  from 

the  solution.  (Equation  (6)  is  adopted  from  [17].) 

The  post-solution  residuals  can  also  be  considered 
to  be  vector  quantities: 

ft  -  djX(Pj  -Xi)xdj  (7) 

with  the  vector  ft  extending  from  the  point  X;  to 
the  nearest  point  on  its  associated  LOP.  The  length 
of  the  vector  is  ft  (since  I  d;  1  =  1).  Considered  as  a 
time  series,  the  vectors  ft.  contain  all  the  informa¬ 
tion  on  the  influence  of  the  observations  on  the  so¬ 
lution,  with  the  effect  of  each  observation  propor¬ 
tional  to  <5?  (or,  if  weighted,  Aft)2). 

Prior  to  a  solution  it  is  easy  to  obtain  an  estimate 
of  its  accuracy  by  considering,  in  general  terms,  the 
geometry  of  the  LOPs.  For  this  purpose,  we  assume 
that  all  the  observations  are  of  similar  quality  and 
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are  well  distributed  in  direction  and  time,  that  is, 
that  there  is  little  or  no  geometric  dilution  of  preci¬ 
sion.  As  shown  in  Figure  2,  each  LOP  is  defined  by 
both  its  “anchor”  in  3-D  space,  a  point  at  the 
assumed  coordinates  of  the  object  observed,  and  its 
direction,  defined  by  the  observation  itself.  Given 
that  the  equation  of  the  LOP  is  X  =  P  +  rd  (where 
r  is  a  scalar  of  arbitrary  value),  the  statistical 
uncertainties  at  a  distance  r  from  the  object  are 
related  by 

Ox  —  Gp  +  r2  (8) 

where  each  a  is  the  root-sum-square  of  the  uncer¬ 
tainties  in  the  respective  vector  components.  Since 
d  is  always  a  unit  vector,  oa  represents  an  angular 
uncertainty  in  radians,  which  is  developed  below; 
we  anticipate  that  it  will  be  closely  related  to  the 
centroiding  error  of  the  imaging  system.  Figure  3 
represents  a  geometric  interpretation  of  (8).  The 
first  term  on  the  right  of  the  equation  represents 
the  average  radius  of  an  ellipsoid  of  uncertainty 
around  the  assumed  position  of  the  observed  object 
due  to  the  likely  errors  in  its  coordinates.  The 
r-term  represents  a  cone  of  expanding  uncertainty 
with  its  axis  along  d,  its  apex  at  the  assumed  posi¬ 
tion  of  the  object  (where  r=0),  and  its  apex  angle 
equal  to  2od.  The  LOP  could  therefore  plausibly  be 
any  line  originating  within  the  ellipsoid  of  uncer¬ 
tainty  with  a  direction  parallel  to  any  line  within 
the  cone  of  uncertainty. 

We  expect  the  observer  to  be  somewhere  within, 
or  not  far  outside  of,  each  LOP’s  volume  of  uncer¬ 
tainty.  So  for  observation  i,  we  have  <5;  «  axirf),  for 
which  we  need  at  least  a  crude  estimate  of  r,,  the 


distance  of  the  object  from  the  observer.  If  the  val¬ 
ues  of  the  <5;’s  computed  this  way  are  similar,  which 
would  generally  be  the  case  if  the  distances  to  the 
observed  objects  were  not  too  different,  then  a  typi¬ 
cal  say  <5,  will  be  a  predictor  of  the  scatter  in  the 
post-fit  residuals.  That  is,  fi2  should  approximately 
equal  the  variance  of  the  fit.  We  therefore  have  a 
simple  way  to  anticipate  the  accuracy  that  can  be 
obtained  by  various  observing  schemes. 

The  one  remaining  piece  of  unfinished  business  is 
determining  the  value  of  ad  to  use  in  Equation  (8). 
If  the  measurement  of  angles  were  absolute,  that 
is,  obtained  from  the  pointing  of  the  imaging  sys¬ 
tem  (say,  from  shaft  encoders  on  the  axes),  then  od 
would  simply  be  the  larger  of  the  mechanical  point¬ 
ing  resolution  or  the  image  resolution.  However, 
such  a  system  would  require  a  transformation  of 
the  angular  measurements  from  an  instrumental 
system  to  a  geodetic  system  before  Equation  (3) 
could  be  applied.  Generally,  that  transformation  is 
unknown.  The  advantage  of  differential  measure¬ 
ments  (i.e.,  nearby  points  measured  with  respect  to 
more  distant  ones  within  a  limited  field  of  view)  is 
that  they  avoid  the  problem  of  the  unknown  instru¬ 
mental  attitude.  The  orientation  of  the  instrument 
need  not  even  be  measured. 

For  differential  measurements,  the  geometry  of 
how  the  observation  vector  is  formed  is  shown  in 
Figure  4.  The  figure  shows  the  effect  of  finite  imag¬ 
ing  resolution,  which  creates  an  ambiguity  in  which 
vector  to  choose.  That  ambiguity  defines  the  angu¬ 
lar  uncertainty  of  the  observation,  od,  used  in 
Equation  (8).  The  figure  is  based  on  the  assumption 
that  that  the  finite  imaging  resolution  is  the  pri¬ 
mary  source  of  uncertainty  and  that  the  errors  in 
the  object  coordinates  are  not  the  limiting  factor. 
The  geometry  yields  the  result  that  ad  — 
A0(| +  ?£?)>  w^ere  A 8  is  the  imaging  resolution,  and 
r  and  r'  are  the  distances  to  the  near  and  far 
objects,  respectively  (we  assume  that  the  far  objects 
have  some  typical,  or  at  least  minimum,  distance). 
We  see  that  as  r'  ->  oo,  od  ->  A0/2,  so  that  very  dis¬ 
tant  background  objects  are  preferred.  Note  that  in 
any  case  we  require  that  the  background  objects  be 
significantly  behind  those  in  the  foreground  so  that 
the  term  does  not  blow  up.  If  a  moving  observer 
waits  for  a  pre-selected  foreground  object  to  appear 
to  line  up  with  a  pre-selected  background  object, 
the  situation  is  similar.  The  diagram  for  this  case  is 
somewhat  different  but  the  resulting  expression  for 
Od  is  the  same. 

EFFECTS  OF  DEVIATIONS  FROM  THE 
IDEAL  TRACK 

Equation  (3)  is  based  on  an  observer  trajectory 
that  curves  only  toward  the  center  of  the  Earth, 
that  is,  the  path  over  ground  is  a  great  circle. 
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Fig.  4-Possible  observation  vectors  for  foreground  objects  viewed  against  background  objects,  at 
distances  from  the  observer  of  r  and  r1,  respectively.  (Not  all  possible  vectors  are  shown.)  The  imag¬ 
ing  system  cannot  distinguish  among  objects  in  the  gray  cone,  defined  by  the  image  resolution  A0. 
The  figure  shows  a  spread  of  possible  observation  vector  directions  (=  LOP  directions)  within  an 
angle  2od. 


Clearly,  this  ideal  will  seldom  play  out  in  practice,  For  example,  it  is  worthwhile  investigating  how 

and  it  is  important  to  evaluate  the  consequences  of  the  algorithm  responds  to  low-grade  accelerations 

deviations  from  the  modeled  motion.  that  might  result  from  a  systematic  change  in  wind 

First,  both  ships  and  aircraft  often  follow  rhumb  or  current  during  the  time  period  covered  by  the 

lines  rather  than  great  circle  routes.  Over  long  dis-  observations.  For  these  cases,  the  results  of  [18] 

tances,  a  great  circle  route  may  be  approximated  provide  some  insight.  The  appendix  of  that  paper 

by  a  series  of  rhumb  lines.  A  rhumb  line  (loxo-  develops  the  general  case  of  a  least-squares  fit  of  a 

drome;  a  track  of  constant  azimuth)  is  a  straight  linear  motion  model  to  position  observations  when 

line  on  a  Mercator  map,  but  it  generally  has  curva-  the  actual  motion  involves  a  weak  acceleration.  The 

tures  in  two  directions  when  viewed  in  a  3-D  coordi-  paper  provides  expressions  for  the  systematic 

nate  system.  Curvature  in  the  horizontal  plane  has  errors  of  the  solved-for  parameters  as  well  as  for 

not  been  accounted  for  in  the  development  here.  the  statistics  of  the  residuals.  Here,  the  residuals 

Rhumb  lines  diverge  most  rapidly  from  great  circles  we  are  interested  in  are  the  distances  between  the 

for  east-west  tracks  (except  for  latitudes  within  a  solution  track  and  the  true  (accelerated)  track, 

few  degrees  of  the  equator,  where  the  divergence  is  evaluated  at  the  times  of  the  observations.  Of 

small  in  any  case)  with  the  effect  being  greater  at  course,  these  residuals  are  only  available  for  nu- 

higher  latitudes.  At  40  deg  latitude,  the  maximum  merical  test  cases,  and  their  statistical  properties 

horizontal  difference  between  a  rhumb  line  and  a  may  differ  from  the  observation  residuals  described 

great  circle  is  41  m  over  a  50  km  track  and  164  m  by  Equation  (6).  Although  the  type  of  observations 

over  a  100  km  track,  if  the  end  points  are  the  same.  and  the  form  of  the  solution  is  different  in  [18]  from 

Except  for  possible  applications  involving  high-  that  considered  here,  numerical  tests  verified  that 

speed  aircraft,  then,  the  systematic  error  of  follow-  the  geometry  is  similar  enough  to  provide  usable 

ing  a  constant  heading  rather  than  a  great-circle  estimates  of  the  magnitude  of  the  effects.  These 

route  would  not  appear  to  be  a  major  issue.  Other  tests  indicated  that  Equation  (3)  generally  produces 

systematic  shifts  in  the  vehicle’s  track  are  likely  to  smaller  increases  in  the  maximum  residual  than 

be  more  important.  predicted  by  [18]  but  somewhat  greater  increases  in 
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the  RMS  of  the  residuals,  relative  to  a  solution 
without  acceleration  included  in  the  vehicle  track. 

For  reasonable  cases  the  solution  adjusts  itself 
such  that  the  residuals  do  not  increase  dramatically 
when  an  acceleration  is  introduced.  The  aforemen¬ 
tioned  paper  predicts  that  the  maximum  residual 
will  increase  by  about  aA/2/12,  where  a  is  the  accel¬ 
eration  and  A t  —  tn  -  t1  is  the  span  of  observation 
times.  For  example,  a  ship  traveling  at  50  km/h  (27 
kn)  will  cover  8.3  Ion  in  10  minutes;  a  100  km/h2 
acceleration  will  shift  the  track  by  1.4  km  over  the 
same  time.  (This  acceleration,  approximately  1  X 
10-3  g,  is  equivalent  to  moving  into  a  9  kn  current 
during  that  time.)  Numerical  simulations  of  this 
scenario  (among  others)  were  computed,  both  with 
and  without  acceleration,  all  with  great-circle  track 
solutions  formed  according  to  Equation  (3)  based  on 
ten  noisy  artificial  data  points.  In  one  typical  case, 
the  maximum  residual  increased  by  0.20  km  (from 
0.29  to  0.49  km)  when  the  acceleration  was  intro¬ 
duced.  The  above  expression  from  [18]  predicted  an 
increase  of  0.23  km.  The  RMS  of  the  residuals 
increased  by  36%  (from  0.19  to  0.26  km);  the  pre¬ 
dicted  increase  was  14%. 

Although  the  past  positions  computed  from  the 
accelerated  solution  were  only  moderately  degraded, 
relative  to  those  from  an  unaccelerated  solution, 
the  solved-for  velocity  does  not  provide  an  accurate 
prediction  of  the  future  track,  whether  it  is  acceler¬ 
ated  or  not.  For  the  above  case,  with  acceleration 
included,  the  difference  between  the  solved-for  ve¬ 
locity  (essentially,  the  average  velocity  over  the  10- 
minute  period  of  the  observations)  and  the  instanta¬ 
neous  velocity  at  either  endpoint  was  8.1  km/h 
(the  prediction  was  8.3  km/h),  that  is,  a  17%  error, 
with  the  error  vector  parallel  to  the  direction  of  the 
acceleration.  Overall,  if  the  solution  were  used  to 
extrapolate  the  observer’s  position  10  minutes 
beyond  the  span  of  observations,  we  should  expect  a 
systematic  error  of  about  1.6  km  if  the  future  motion 
was  unaccelerated  and  3.0  km  if  the  acceleration 
continued.  On  the  other  hand,  if  both  the  accelera¬ 
tion  and  the  observations  continue,  a  set  of  rolling 
solutions  (i.e.,  using  observations  within  a  moving 
window  of  time)  could  provide  the  acceleration  from 
the  continuous  change  in  velocity  from  one  solution 
to  the  next. 

We  can  construct  approximate  formulas  for  the 
applicability  of  the  algorithm  in  the  presence  of 
acceleration,  based  on  either  of  two  criteria.  We  use 
the  expressions  in  [18]  that  show  that  at  the  begin¬ 
ning  and  end  of  the  observation  interval,  which 
spans  time  At,  the  solution’s  velocity  will  differ 
from  the  observer’s  instantaneous  velocity  by  a  At 1 2 
and  the  systematic  error  in  position  there  is  a  At2! 
12.  The  uncertainty  of  the  solution  projected  to 

time  t  is  o{t)  —  +  ( av(t  —  to))2,  where  ax  and  av 


are  the  formal  uncertainties  in  position  and  veloc¬ 
ity,  respectively  (their  covariance  is  ignored  here), 
and  t0  is  the  epoch  chosen  for  the  solution  parame¬ 
ters,  which  we  assume  to  be  at  the  end  of  the  obser-  | 

vations.  Then  the  expected  position  errors,  e,  at 
times  t0  and  t0  +  At  are,  respectively  f 


e(f0)  =^~aAt2  +  ax  and 

(9) 

e(to  +  At)  —  Y2a&t2  +  \j Gx  +  i^vAt)2 

In  the  latter  case,  we  have  assumed  that  the  accel¬ 
eration  does  not  continue  beyond  the  end  of  the 
observations;  if  the  acceleration  continues,  the  frac¬ 
tion  yi  would  be  replaced  by  ||  If  emax  is  the  maxi¬ 
mum  allowable  navigation  error,  then  using  the 
above  expressions,  we  find  that  as  long  as 


the  solution  will  be  acceptable.  “Acceptable”  means, 
in  the  first  case,  that  the  expected  position  error  at 
the  end  of  the  observation  interval  (at  /o)  will  be 
less  than  emax,  and  in  the  second  case  that  the  solu¬ 
tion  could  be  projected  another  interval  At  into  the 
future  (to  to  +  At)  with  error  less  than  emax.  In  the 
latter  expression,  m  —  12/13  if  the  acceleration  con¬ 
tinues  and  m  —  12/7  if  it  does  not.  In  the  example 
referred  to  in  the  preceding  paragraphs,  ax  =  0.25 
km  and  av  =  2.24  km/h  for  the  solution  with  accel¬ 
eration,  and  At  —  1/6  h.  If  we  assume  that  the  accel¬ 
eration  continues  after  the  observations  (t  >  to),  so 
that  m  =  12/13,  and  we  set  emax  =  1  km,  then  the 
acceleration  a  would  have  to  be  less  than  18  km/h2 
(equivalent  to  moving  into  a  1.7  kn  current  over  10 
minutes)  for  the  solution  to  be  acceptable  based  on 
its  predictive  ability  beyond  the  observations.  But  if 
only  the  positional  accuracy  at  the  end  of  the  obser¬ 
vation  span  were  the  criterion,  the  acceleration 
could  be  as  high  as  400  km/h2,  i.e.,  over  20  times 
greater,  for  the  solution  to  be  acceptable. 

Of  course,  the  real  world  is  more  complicated.  In 
most  cases  we  would  have  little  information  on  the 
magnitude  of  any  such  acceleration,  and  it  is 
unlikely  to  be  constant  for  extended  periods  of  time 
anyway  The  actual  track  of  the  vehicle  will  in  gen¬ 
eral  be  subject  to  both  systematic  and  random-walk 
shifts  due  to  changing  wind  or  currents  or  steering 
or  propulsion  variations.  The  navigation  solution 
described  here,  although  computed  for  the  specific 
instant  to,  really  reflects  a  kind  of  average  track  of 
the  vehicle  during  the  observations.  In  that  way  it 
differs  from  near-instantaneous  determinations  of 
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position  and  velocity  such  as  from  GPS.  The  accu¬ 
racy  of  dead-reckoning  predictions  based  on  a  single 
navigation  solution  of  either  kind  will  vary  widely 
depending  on  conditions.  Although  the  Equation  (3) 
algorithm  has  limited  predictive  ability  when  a  con¬ 
stant  acceleration  is  present,  it  would  do  well  if  the 
track  variations  were  stochastic  or  nearly  so  and 
enough  observations  were  used  to  provide  an 
unbiased  sample.  Stochastic  track  variations  have 
the  same  general  effect  on  a  solution  as  an  increase 
in  the  random  errors  of  observation. 


THE  USE  OF  SATELLITE  OBSERVATIONS 

As  was  mentioned  in  the  second  section  of  this 
paper,  a  specific  proposed  application  would  use  op¬ 
tical  or  near-infrared  measurements  of  the  angular 
positions  of  satellites  observed  against  a  star  back¬ 
ground.  Such  a  system  could  provide  a  standalone 
backup  against  jamming  of  GPS  signals,  although 
likely  of  lesser  accuracy.  GPS  satellites  present 
point-like  images  (about  0.1  arcsecond  across,  less 
than  the  atmospheric  “seeing”  disk)  in  the  visual 
magnitude  range  of  11-14,  depending  on  the  Sun- 
satellite-observer  angle  [19].  The  main  challenge  in 
obtaining  such  observations  is  the  large  difference 
in  the  angular  motions  of  the  GPS  satellites  and 
the  background  stars,  which  can  exceed  30  arcsec/s 
(0.5  deg/min)  as  seen  from  the  surface  of  the  Earth. 
Nevertheless,  the  satellites  are  observable  in  small 
telescopes  with  modern  electronic  sensors.  Because 
the  observed  satellites  are  at  finite  distances,  with 
geocentric  coordinates  known  to  a  meter  or  better 
(readily  available  on  the  Internet)  a  straightfor¬ 
ward  triangulation  method,  such  as  the  one  pre¬ 
sented  in  this  paper,  is  feasible.  Other  satellites 
with  accurately  known  orbits,  e.g.,  geosynchronous 
communications  satellites,  or  low-Earth-orbit  (LEO) 
geodetic  satellites,  might  also  serve  as  observational 
targets. 

Of  course,  the  observer  must  know  his  position  at 
least  approximately  to  be  able  to  plan  observations. 
For  GPS  satellites,  which  are  at  a  minimum  height 
of  about  20,000  km,  the  observer’s  position  must  be 
known  only  to  within  175  km  for  a  satellite  to  fall 
within  a  pre-pointed  1  deg  field  of  view  (angular 
error  0.5  deg).  For  LEO  geodetic  satellites,  with 
typical  heights  of  1000  km,  the  observer’s  a  priori 
position  must  be  known  to  within  9  km.  For  cam¬ 
eras  with  apertures  of  about  20  cm,  and  charge- 
coupled  device  (CCD)  sensors  in  the  focal  plane,  ex¬ 
posure  times  would  be  less  than  a  minute,  often 
only  a  few  seconds,  depending  on  the  satellite  ge¬ 
ometry.  (The  required  exposure  is  thus  less  than 
the  time  needed  for  the  satellite  to  cross  a  1  deg 
star  field.)  The  CCDs  could  work  in  either  conven¬ 
tional  “stare”  mode  or  time-delay  integration  (TDI) 


mode.  In  TDI  mode,  charge  is  moved  across  the  sen¬ 
sor  chip  at  a  controlled  rate  that  matches  the  angu¬ 
lar  motion  of  the  object  of  interest,  so  that  the 
object’s  image  accumulates  charge  as  it  moves 
across  the  chip.  A  camera  designed  to  observe  GPS 
satellites  would  probably  use  both  types  of  CCDs  in 
the  focal  plane,  one  type  for  the  satellite,  the  other 
for  the  stars.  Image  timing  accurate  to  about  a 
millisecond  in  UTC  is  required.  Generally,  GPS  sat¬ 
ellites  would  be  night  objects,  but  some  daylight 
observations  might  be  possible  with  near-infrared 
sensors.  Successful  daytime  experimental  observa¬ 
tions  of  LEO  satellites  have  been  carried  out  with 
such  sensors. 

Each  satellite  observation  would  be  expressed  in 
some  sort  of  “space  fixed”  celestial  coordinate  sys¬ 
tem  while  the  satellite  ephemeris  position  would 
most  likely  be  expressed  in  a  geodetic  system  that 
rotates  with  the  Earth;  a  transformation  of  the 
observational  data  from  the  celestial  to  the  terres¬ 
trial  system  will  be  required.  In  practice,  there  are 
two  fundamental  coordinate  systems  that  are  most 
appropriate  for  this  problem:  the  International  Ter¬ 
restrial  Reference  System  (ITRS)  and  the  Geocen¬ 
tric  Celestial  Reference  System  (GCRS).  The  ITRS 
is  a  3-D  geocentric  system  that  rotates  with  the 
Earth,  so,  loosely  speaking,  it  is  a  “crust-fixed”  sys¬ 
tem;  more  precisely,  its  axes  have  no  net  rotation 
with  respect  to  the  ensemble  of  ITRS  defining  sta¬ 
tions.  Geodetic  positions  in  the  ITRS  coincide  (within 
several  cm)  with  positions  measured  with  respect  to 
the  WGS-84  ellipsoid,  e.g.,  from  GPS.  The  ITRS  is 
also  referred  to  in  GPS  literature  as  the  Earth-cen¬ 
tered  Earth-fixed  (ECEF)  system.  The  GCRS  is  the 
geocentric  equivalent  of  the  International  Celestial 
Reference  System  (ICRS),  which  has  its  origin  at 
the  solar  system  barycenter  and  which  is  used  for 
modern  star  catalogs  and  planetary  ephemerides. 
The  GCRS  axes  have  no  net  rotation  with  respect  to 
distant  objects  in  the  universe.  The  GCRS  is  a  natu¬ 
ral  system  for  expressing  the  positions  of  stars  as 
they  would  be  seen  by  an  observer  on  or  near  the 
Earth  at  a  specific  time. 

Since  both  the  ITRS  and  GCRS  are  geocentric 
systems,  the  transformation  between  them  consists 
of  a  series  of  rotation  matrices: 

i’gcrs  =  BPNSW  riTRs  (11) 

where  riTRs  is  a  vector  in  the  terrestrial  system 
and  i’gcrs  is  the  corresponding  vector  in  the  celes¬ 
tial  system.  The  matrices  listed  account  for,  from 
right  to  left,  polar  motion,  sidereal  time,  nutation, 
precession,  and  a  constant  “frame  bias.”  These  coor¬ 
dinate  systems  and  the  transformations  that  link 
them  are  described  more  fully  in  [20]  and  there  are 
several  software  packages  available  on  the  web  for 
carrying  out  this  transformation  or  its  inverse.  This 
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transformation  forms  the  link  between  each  satel¬ 
lite  observation  (right  ascension  and  declination 
with  respect  to  the  GCRS)  and  the  ephemeris  posi¬ 
tion  of  the  satellite  at  the  time  of  the  observation 
(X,Y,Z  with  respect  to  the  ITRS).  Actually,  the 
ephemeris  position  of  the  satellite  should  be  ob¬ 
tained  for  the  time  of  observation  minus  the  light¬ 
time  to  the  GPS  satellite  («0.07  s).  The  light-time 
calculation  can  be  fairly  crude  and  could  be  based 
on  the  distance  of  the  satellite  from  the  geocenter 
(rs26,600  km)  and  its  observed  zenith  angle;  alter¬ 
natively,  an  approximate  observer  location  and  the 
satellite  ephemeris  could  be  used.  (An  error  of  750 
km  in  the  light  travel  distance  results  in  only  a  10 
m  error  in  the  satellite’s  position.) 

If  the  image  centroiding  errors  are,  say,  1-2  arc- 
sec,  then,  using  rJ(i  =  Af)/2  =  5  /irad,  uP  —  1m,  and 
r  «  23,000  km,  the  expected  scatter  in  the  triangula¬ 
tion  residuals,  using  Equation  (8),  would  be  around 
115  m.  However,  tracking  a  satellite  over  just  a  few 
minutes  could  provide  a  large  number  of  independ¬ 
ent  observations,  and  the  uncertainty  of  the  position 
solution  could  conceivably  be  reduced  to  a  few  tens 
of  meters.  An  advantage  of  using  GPS  satellites  is 
that  the  constellation  is  designed  to  provide  good 
ranging  geometry  (i.e.,  to  minimize  GDOP),  so  usu¬ 
ally  the  observational  geometry  would  also  be  favor¬ 
able  for  the  kind  of  triangulation  described  in  this 
paper. 

Unlike  traditional  celestial  navigation  where  the 
observed  stars  are  assumed  to  be  infinitely  distant, 
and  triangulation  is  therefore  not  possible,  the  sat¬ 
ellite  scheme  does  not  require  any  reference  to  the 
local  vertical.  This  is  important  for  moving  observ¬ 
ers,  where  a  precise  determination  of  the  local  ver¬ 
tical  cannot  be  made  using  direct  onboard  (“lab”) 
measurements,  since  the  gravity  and  acceleration 
vectors  are  inseparable.  (Inertial  navigation  sys¬ 
tems  can  provide  a  computed  estimate  of  the  local 
vertical.)  The  horizon,  which  ideally  defines  an 
external  plane  orthogonal  to  the  local  vertical,  is  of¬ 
ten  not  visible  or  accurately  measurable.  Even 
when  clearly  visible  and  sharply  defined,  the  hori¬ 
zon  is  actually  a  warped  circle,  at  the  level  of  preci¬ 
sion  needed,  due  to  direction-dependent  low-level 
atmospheric  refraction.  The  satellite  scheme  cir¬ 
cumvents  the  local  vertical  or  horizon  problem; 
additionally,  because  it  uses  small-field  measure¬ 
ments  from  an  electronic  image,  there  is  no  need 
for  precise  large-angle  calibration.  Thus,  for  moving 
observers  that  require  a  supplement  or  backup  to 
ordinary  radio-based  GPS  navigation,  angular 
measurements  of  GPS  satellites  with  respect  to  the 
stars  have  several  fundamental  advantages  over 
traditional  celestial  navigation,  even  if  the  latter  is 
automated. 

Additionally,  the  satellites,  like  the  stars,  provide 
an  absolute  attitude  reference.  If  the  observations 


used  for  the  position  and  velocity  solution  can  be 
expressed  in  instrumental  coordinates,  then  we 
have  two  bases  for  each  such  vector:  an  external 
reference  system  (either  the  ITRS  or  the  GCRS) 
and  the  instrumental  reference  system.  That  is  the 
information  needed  to  solve  “Wahba’s  problem” 
[21],  about  which  there  is  an  extensive  literature 
and  software  base,  and  determine  the  attitude  of 
the  instrument. 


NUMERICAL  SIMULATIONS 

A  software  package  was  created  to  test  the  math¬ 
ematics  presented  in  this  paper  and  to  explore  the 
properties  of  the  solutions.  The  software  computes 
the  track  of  a  hypothetical  vehicle  across  the  sur¬ 
face  of  the  Earth,  given  an  initial  date  and  time 
and  values  for  latitude,  longitude,  course,  and 
speed;  either  a  great-circle  or  a  rhumb-line  track 
can  be  selected.  This  track  is  “truth.”  The  user  can 
specify  a  time  span  within  which  a  selected  number 
of  artificial  observations  will  be  created.  For  each 
observation,  the  software  identifies  a  target  object, 
obtains  its  coordinates,  and  computes  an  observa¬ 
tion  vector,  with  the  target  coordinates  and  “ob¬ 
served”  direction  subject  to  random  errors.  The  en¬ 
semble  of  observation  times,  target  coordinates, 
and  observation  vectors  is  then  sent  to  the  routine 
that  sets  up  and  solves  Equation  (3).  The  solution 
yields  the  vehicle’s  position  (at  a  pre-selected  time) 
and  velocity  in  geocentric  rectangular  coordinates, 
which  can  be  transformed  into  latitude,  longitude, 
height,  course,  speed,  and  rate  of  change  of  height. 
This  allows  the  computation  of  a  “solution  track,” 
which  is  compared  to  the  “true  track”  for  the  span 
of  time  covered  by  the  observations;  the  differences 
are  sampled  at  the  times  of  the  observations.  For 
simplicity,  a  spherical  Earth  with  a  radius  of 
6378.137  1cm  was  assumed  for  all  the  tests  reported 
here.  A  computer  script  was  written  that  allowed 
large  numbers  of  solutions  to  be  formed  and  com¬ 
pared  for  a  given  test  track,  using  a  Monte-Carlo- 
type  scheme.  Many  different  solution  scenarios 
were  tested  this  way. 

In  computing  the  artificial  observations,  the  soft¬ 
ware  attempts  to  spread  the  observation  times 
more  or  less  evenly  across  the  specified  time  span 
(say,  a  half-hour),  although  the  exact  times  chosen 
are  random.  For  each  time  chosen,  the  software 
checks  the  elevation  angle  and  azimuth  of  GPS  sat¬ 
ellites  as  they  would  appear  from  the  vehicle  at 
that  time.  Actual  GPS  positions  are  used,  derived 
from  NORAD  two-line  orbital  elements  propagated 
to  the  time  of  observation  using  SPACETRACK 
software  [22],  Any  satellite  to  be  “observed”  must 
be  at  least  30  deg  above  the  horizon,  and  the  soft¬ 
ware  attempts  to  select  a  satellite  that  significantly 
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Table  1  — Hypothetical  True  Track  and  Artficial  Observations 


UTC  Time 

Distance 

km 

Latitude 

deg 

Longitude 

deg 

Obs 

# 

Sat  ID 

Altitude 

deg 

Azimuth 

deg 

Range 

km 

03:31:21 

-23.9 

+39.89 

-50.24 

1 

PRN  06 

32.3 

146.1 

22635 

03:39:27 

-17.1 

+39.92 

-50.17 

2 

PRN  01 

84.9 

264.9 

20315 

03:45:57 

-11.7 

+39.95 

-50.12 

3 

PRN  30 

40.8 

49.6 

21678 

03:47:43 

-10.2 

+39.95 

-50.10 

4 

PRN  14 

54.1 

150.6 

21106 

03:51:05 

-7.4 

+39.97 

-50.08 

5 

PRN  01 

83.1 

208.0 

20315 

03:52:37 

-6.2 

+39.97 

-50.06 

6 

PRN  31 

67.0 

334.5 

20742 

03:57:09 

-2.4 

+39.99 

-50.02 

7 

PRN  29 

42.0 

90.0 

21913 

03:59:15 

04:00:00 

-0.6 

0.0 

+40.00 

+40.00 

-50.01 

-50.00 

8 

PRN  01 

79.6 

191.8 

20348 

differs  in  azimuth  from  those  recently  observed.  (In 
a  real-world  application,  satellite  positions  com¬ 
puted  by  SPACETRACK,  based  on  two-line  orbital 
elements,  would  not  be  nearly  accurate  enough  for 
this  purpose.  GPS  positions  with  errors  of  a  meter 
or  less  would  have  to  be  obtained  from  Interna¬ 
tional  GPS  Service  data  sets  distributed  daily  on 
the  IGS  web  site.) 

Each  observation  vector  —  that  is,  the  unit  vector 
that  defines  the  direction  of  the  LOP  —  is  computed 
from  the  difference  between  the  GPS  position  and 
the  vehicle  position  for  the  same  instant,  both 
expressed  in  an  ITRF-like  rectangular  coordinate 
system  that  rotates  with  the  Earth  (i.e.,  an  ECEF 
system).  No  attempt  is  made  to  simulate  an  actual 
observing  system,  such  as  a  telescopic  CCD  camera, 
or  to  determine  which  stars  the  satellite  would  be 
seen  against.  Light-time  is  neglected,  as  are  effects 
such  as  refraction  and  aberration  that  would  be  the 
same  for  the  satellite  and  the  background  stars. 
However,  each  observation  is  modified  by  adding  a 
random  error  in  angle,  in  a  random  direction.  The 
magnitudes  of  the  angular  errors  are  normally  dis¬ 
tributed  with  a  standard  deviation  (in  arcseconds) 
provided  by  the  user.  Weights  are  not  assigned  to 
the  observations. 

The  satellite  positions  are  also  subject  to  random 
position  errors.  Just  as  for  the  observation  errors, 
the  user  provides  a  standard  deviation  (in  kilo¬ 
meters)  for  the  satellite  position  errors,  which  are 
then  normally  distributed  in  each  of  three  dimen¬ 
sions.  (Real  satellite  position  errors  are  generally 
higher  in  the  along-track  than  cross-track  direc¬ 
tions.) 

A  few  generalities  about  the  solutions:  When  the 
user-selected  obsei'vational  errors  (standard  devia¬ 
tions)  are  set  to  zero,  the  solutions  are  for  practical 
purposes  exact  for  great-circle  tracks  of  a  few  tens 
of  kilometers  long,  as  expected.  That  is,  the  lati¬ 
tude,  longitude,  course,  and  speed  obtained  from 
the  solution’s  position  and  velocity  vectors  are 
essentially  identical  to  those  that  describe  the  true 
track  for  the  same  instant.  Along  the  solution  track, 
the  positional  errors  are  less  than  1  m  compared  to 


the  true  track.  As  the  time  span,  hence  the  dis¬ 
tance,  over  which  observations  are  obtained 
increases,  the  errors  for  such  “perfect  observation” 
test  cases  also  gradually  increase,  due  to  the  small 
errors  in  the  curvature  term  described  in  [14].  Most 
of  these  solutions  used  the  correct  vehicle  speed  v 
in  the  curvature  term.  However,  even  if  v  has  a  rel¬ 
atively  large  error  (say,  30%),  it  is  only  necessary  to 
solve  Equation  (3)  twice,  using,  the  second  time, 
the  value  of  v  obtained  from  the  velocity  vector 
from  the  first  solution.  Tests  showed  that  the  sec¬ 
ond  solution  is  essentially  identical  to  a  single  solu¬ 
tion  with  v  known  exactly. 

Example:  A  typical  solution  for  a  ship’s  rhumb¬ 
line  track  is  shown  in  Tables  1,  2,  and  3.  The  ship 
was  hypothesized  to  be  moving  at  a  constant  speed 
of  50  km/h  (27  kn)  along  a  course  60  deg  from  true 
North.  On  2008  February  19  at  04:00  UTC,  the 
ship’s  position  was  50  deg  West  longitude,  40  deg 
North  latitude.  Eight  artificial  observations  of  vari¬ 
ous  GPS  satellites  were  computed  for  times  distrib¬ 
uted  over  the  previous  half-hour,  corresponding  to  a 
25-km  track  at  the  ship’s  speed.  The  standard  devi¬ 
ation  of  the  observational  errors  was  set  at  1  arcsec 
(5  /trad)  and  the  standard  deviation  of  the  GPS 
coordinates  was  set  at  5  m.  Equation  (8)  predicts  a 
scatter  in  the  residuals  of  just  over  100  m  for  this 
case,  most  of  which  is  from  the  angular  uncertainty 
of  the  observations.  The  solution  was  computed  for 
04:00  UTC.  As  with  all  solutions  to  (3),  only  three 
data  elements  were  involved  for  each  observation: 
the  time,  the  observed  satellite’s  direction  vector, 
and  the  observed  satellite’s  geocentric  position  vec¬ 
tor.  Both  vectors  were  expressed  in  ITRF-like  rec¬ 
tangular  coordinates  that  rotate  with  the  Earth.  (In 
Table  1,  the  range  to  each  satellite  is  shown  for  in¬ 
formation  only  and  was  not  used  in  the  solution.) 
Note  that  three  observations  of  PRN  01  were 
selected  as  it  passed  nearly  overhead;  although  the 
observing  geometry  is  similar  in  each  case,  they 
sample  the  ship’s  position  at  different  times. 

The  solution  parameters  are  given  in  Table  2. 
The  top  half  of  the  table  gives  the  position  and  ve¬ 
locity  vectors  expressed  in  geocentric  rectangular 
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Table  2— Solution  for  Ship’s  Position  and  Velocity  at  04:00  UTC 


X 

km 

y 

km 

Z 

km 

A 

km/h 

I 

km/h 

km/h 

Solution 

4-  3140.644 

-3742.714 

+4099.755 

+22.879 

+40.496 

+19.031 

Truth 

+  3140.619 

-3742.844 

+4099.787 

+22.841 

+40.144 

+19.151 

Difference 

+0.025 

+0.131 

-0.032 

+0.038 

+0.352 

-0.120 

Formal  Mean  Error 

0.098 

0.078 

0.082 

0.436 

0.311 

0.273 

Latitude 

Longitude 

Height 

Course 

Speed 

h 

deg 

deg 

km 

deg 

km/h 

km/h 

Solution 

+40.00026 

-49.99879 

-0.085 

60.0814 

50.254 

-0.264 

Truth 

+40.00000 

-50.00000 

0.000 

60.0000 

50.000 

0.000 

Difference 

+0.00026 

+0.00121 

-0.085 

+0.0814 

+0.254 

-0.264 

Table  3 — Post-Solution  Residuals 


Obs  # 

Distance 

km 

dRy 

km 

dNr 

km 

dEy 

km 

dUr 

km 

dRi 

1cm 

dNz, 

1cm 

dEL 

km 

dUx. 

km 

1 

-23.9 

0.051 

-0.030 

-0.003 

+0.041 

0.044 

+0.025 

-0.004 

+0.036 

2 

-17.1 

0.025 

-0.008 

+0.024 

+0.006 

0.132 

-0.118 

-0.057 

-0.006 

3 

-11.7 

0.053 

+0.007 

+0.047 

-0.023 

0.101 

+0.088 

-0.029 

-0.041 

4 

-10.2 

0.063 

+0.011 

+0.054 

-0.031 

0.182 

+0.047 

+0.172 

-0.031 

5 

-7.4 

0.082 

+0.017 

+0.066 

-0.046 

0.049 

-0.017 

+0.046 

+0.001 

6 

-6.2 

0.092 

+0.020 

+0.073 

-0.052 

0.093 

+0.006 

-0.091 

-0.019 

7 

-2.4 

0.119 

+0.026 

+0.091 

-0.072 

0.097 

+0.025 

-0.063 

+0.069 

8 

RMS 

-0.6 

0.132 

0.090 

+0.029 

+0.100 

-0.082 

0.061 

0.111 

-0.055 

+0.024 

-0.009 

: 

| 


\ 


coordinates,  directly  from  the  solution.  The  lower 
half  of  the  table  converts  that  data  to  geodetic  coor¬ 
dinates  on  the  spherical  Earth  model  used.  The  last 
column  at  lower  right  is  the  rate  of  change  of 
height. 

Table  3  shows  the  post-solution  residuals  for  each 
observation.  The  residuals  on  the  left  side  of  the  ta¬ 
ble  indicate  the  distance  between  the  solution  track 
and  the  true  track  at  the  time  of  each  observation; 
dRr  is  the  total  distance  between  tracks,  while  dNr, 
dET,  and  dUy  are  the  components  of  that  distance  in 
the  topocentric  directions  North,  East,  and  up.  The 
run  of  dRy  values  indicates  the  overall  position 
error  of  the  solution.  The  residuals  on  the  right 
side  of  the  table  are  similar,  but  they  indicate  the 
distance  between  the  solution  track  and  the  line 
of  position  defined  by  each  observation.  As  expected, 
the  line-of-position  residuals,  which  are  the  basis 
for  the  solution,  appear  randomly  distributed 
whereas  the  track  residuals  show  systematic  trends. 
It  is  the  RMS  of  the  line-of-position  residuals  that 
are  predicted  by  (8).  Note  that  the  dUy  track  resid¬ 
uals  show  the  effect  of  small,  non-zero  values  for  the 
height  and  height-rate  parameters;  no  attempt  was 
made  to  constrain  the  solution’s  track  to  the  surface 
of  the  Earth.  The  dNy  and  dEy  residuals  include  the 
rhumb-line  versus  great-circle  errors,  although  they 
are  small  (<10  m)  for  this  case. 


One  hundred  solutions  for  the  same  hypothetical 
rhumb-line  track  were  computed,  each  with  either 
6,  7,  8,  or  9  observations  spread  over  the  same  half- 
hour  (25  solutions  each).  The  median  position  error  : 

(distance  dRy)  was  70  m,  and  74%  of  the  position 
errors  were  less  than  100  m.  The  distribution  of  the 
position  errors  appeared  Poisson-like,  with  only  one 
error  out  of  750  samples  over  300  m  (387  m)  and 
just  25  (3%)  over  200  m.  In  these  tests,  there  was 
only  a  very  slight  hint  that  the  solutions  with  the 
greater  number  of  observations  were  better  overall 
than  those  with  fewer;  generally  the  solution-to- 
solution  variations  dominated  the  statistics.  Other 
numerical  tests  indicated  that  significantly  increas-  } 

ing  the  number  of  observations  does  have  a  measur¬ 
able  effect  in  reducing  the  track  errors.  i 


A  unique,  closed-form  algorithm  has  been  pre¬ 
sented  that  provides  a  3-D  navigational  solution  for 
position  and  velocity  given  a  sequence  of  angles- 
only'  measurements.  The  algorithm  includes  a  cor¬ 
rection  term  for  the  curvature  of  the  Earth,  so  that 
observations  can  be  collected  over  extended  tracks 
if  necessary.  The  kind  of  measurements  proposed 
for  use  with  this  algorithm  are  those  in  which  the 
apparent  positions  of  foreground  objects  are  imaged 


CONCLUSION 


and  measured  against  background  objects,  and 
coordinates  in  some  well-defined  reference  system 
(or  systems)  are  known  for  both.  This  represents  a 
large  class  of  applications,  which  includes  the  ob¬ 
servation  of  satellites  with  accurate  ephemerides, 
such  as  GPS,  observed  against  a  star  background. 

Numerical  simulations  have  shown  the  algorithm 
to  be  mathematically  correct  and  the  solutions  ro¬ 
bust,  even  when  the  number  of  observations  is  small. 
The  solutions  provide  timely  navigation  information 
at  a  useful  level  of  accuracy  for  satellite-observing 
scenarios  that  are  realistic  and  based  on  current 
instrumentation  capabilities  and  data  availability. 

The  method  described  in  this  paper  and  one  ver¬ 
sion  of  a  physical  navigation  system  that  imple¬ 
ments  it  are  the  subject  of  U.S.  patent  applications 
12/552534  and  12/319651  (20090177398),  respec¬ 
tively. 
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Javad  GNSS,  Inc. 

John  Deere 
Kongsberg  Seatex  AS 
KVH  Industries,  Inc. 

L-3  Communications/IEC 
LinQuest  Corporation 
Lockheed  Martin 
Loctronix  Corporation 
Maxtena,  Inc. 

MTRG  CESRE  CSIRO 
National  Instruments 
Navcom  Technology,  Inc. 

Navcon 
NavtechGPS 
Nikon-Trimble  Co.,  Ltd. 
Northrop  Grumman 
Electronic  Systems 
Novariant 
Novatel,  Inc. 

Overlook  Systems  Technologies 

PETROBRAS 

Point,  Inc. 

PREDESA,  LLC 
Rakon  Ltd. 

Raven  Industries-ATC 
Raytheon 

Renesas  Design  France 
Rincon  Research  Corporation 
Rockwell  Collins 
Romona,  Inc. 

RX  Networks  Inc. 

Sandia  National  Laboratories 
Sarantel  Ltd. 

Selex  SI  US 

Sensonor  Technologies  AS 
Septentrio  Satellite  Navigation 
SiRF  Technology,  Inc. 

Skyguide  (Swiss  Air  Navigation 
Services) 

Sokkia  Topcon  Co.,  Ltd.,  Technical 
Administration 
Southern  Avionics  Co. 

Spirent  Communications  SW  Ltd. 


Spirent  Federal  Systems 
ST  Ericsson  UK  Ltd. 

ST  Microelectronics 

Steve  Lieber  &  Associates,  Inc. 

Symmetrieom 

Systran  Donner  Inertial 

Tubitak  Sage 

The  Aerospace  Corporation 
The  MITRE  Corporation 
Topcon  Positioning  Systems 

Trimble 

Trimble  Navigation,  Ltd. 

u-blox  ag 

Universal  Avionics  Corporation 
UrsaNav,  Inc. 

White  Electronic  Designs  Corporation 
Xsens  Technologies  B.V. 

Government  Agencies 

FAA  ATO-OPS  Planning-Systems 
Engrg 

Federal  Aviation  Administration 

National  Geodetic  Survey 

SPAWAR  Systems  Center  Atlantic 

USAF  GPS  Directorate 

USAF  GPS  Wing 

U.S.  Air  Force  746  TS 

U.S.  Coast  Guard  Command  & 

Control  Engineering  Center 
U.S.  Coast  Guard  NAVCEN 
U.S.  Coast  Guard 
USCG  R&D  Center 
USDOT/VNTSC 
U.S.  Naval  Observatory 

Educational  Institutions 

Ohio  University,  Avionics  Engineering 
Center 

PSU/ARL  Navigation  R&D  Center 
The  Johns  Hopkins  University, 
Applied  Physics  Laboratory 
The  University  of  Texas 
U.S.  Coast  Guard  Academy 
U.S.  Merchant  Marine  Academy 
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